function svd_experiment()
% SVD_EXPERIMENT 奇异值分解综合实验
% 
% 本实验涵盖SVD的基本概念、算法和各种应用

fprintf('=== 奇异值分解综合实验 ===\n\n');

%% 实验1：SVD基础理解
fprintf('实验1：SVD基础理解\n');
fprintf('目标：理解SVD的基本概念和几何意义\n\n');

% 简单的2×3矩阵
A1 = [3, 2, 2; 2, 3, -2];
fprintf('矩阵 A1 (%d×%d) =\n', size(A1));
disp(A1);

% 计算SVD
[U1, S1, V1] = svd(A1);
fprintf('U矩阵 (%d×%d):\n', size(U1)); disp(U1);
fprintf('奇异值矩阵 S (%d×%d):\n', size(S1)); disp(S1);
fprintf('V矩阵 (%d×%d):\n', size(V1)); disp(V1);

% 验证分解
reconstruction_error = norm(U1 * S1 * V1' - A1, 'fro');
fprintf('重构误差: %e\n', reconstruction_error);

% 奇异值
singular_values = diag(S1);
fprintf('奇异值: σ₁ = %.6f, σ₂ = %.6f\n', singular_values(1:2));

% 几何解释：单位圆的变换
theta = linspace(0, 2*pi, 100);
unit_circle = [cos(theta); sin(theta)];
ellipse = A1 * [unit_circle; zeros(1, length(theta))];

figure('Name', 'SVD几何解释');
subplot(1, 2, 1);
plot(unit_circle(1, :), unit_circle(2, :), 'b-', 'LineWidth', 2);
axis equal; grid on; title('单位圆');
xlabel('x'); ylabel('y');

subplot(1, 2, 2);
plot(ellipse(1, :), ellipse(2, :), 'r-', 'LineWidth', 2);
axis equal; grid on; title('变换后的椭圆');
xlabel('x'); ylabel('y');

%% 实验2：低秩近似
fprintf('\n实验2：低秩近似\n');
fprintf('目标：理解SVD在低秩近似中的应用\n\n');

% 构造低秩矩阵加噪声
rank_true = 3;
m = 8; n = 6;
U_true = randn(m, rank_true);
V_true = randn(n, rank_true);
A_true = U_true * V_true';
noise = 0.1 * randn(m, n);
A2 = A_true + noise;

fprintf('原始矩阵: %d×%d, 真实秩: %d\n', m, n, rank_true);
fprintf('添加噪声后的矩阵秩: %d\n', rank(A2));

% SVD分解
[U2, S2, V2] = svd(A2);
s2 = diag(S2);
fprintf('奇异值: [%.4f, %.4f, %.4f, %.4f, %.4f, %.4f]\n', s2);

% 不同秩的近似
ranks_to_test = [1, 2, 3, 4, 5];
fprintf('\n低秩近似结果:\n');
fprintf('秩    近似误差    压缩比    相对误差\n');
fprintf('--    --------    ------    --------\n');

for r = ranks_to_test
    % r秩近似
    A_approx = U2(:, 1:r) * S2(1:r, 1:r) * V2(:, 1:r)';
    
    % 计算误差
    approx_error = norm(A2 - A_approx, 'fro');
    relative_error = approx_error / norm(A2, 'fro');
    
    % 压缩比
    original_elements = m * n;
    compressed_elements = r * (m + n + 1);
    compression_ratio = compressed_elements / original_elements;
    
    fprintf('%2d    %8.4f    %6.3f    %8.4f\n', ...
            r, approx_error, compression_ratio, relative_error);
end

%% 实验3：主成分分析(PCA)
fprintf('\n实验3：主成分分析(PCA)\n');
fprintf('目标：使用SVD实现PCA降维\n\n');

% 生成2D数据
n_samples = 100;
% 原始数据在一个方向上有更大的方差
X_original = randn(n_samples, 2);
X_original(:, 2) = 0.3 * X_original(:, 2);  % 压缩y方向的方差

% 旋转数据
angle = pi/4;
R = [cos(angle), -sin(angle); sin(angle), cos(angle)];
X3 = (R * X_original')';

fprintf('数据维度: %d×%d\n', size(X3));
fprintf('数据均值: [%.4f, %.4f]\n', mean(X3));

% 中心化数据
X_centered = X3 - mean(X3);

% 使用SVD进行PCA
[U3, S3, V3] = svd(X_centered, 0);
principal_components = V3;
explained_variance = diag(S3).^2 / (n_samples - 1);

fprintf('主成分方向:\n');
fprintf('PC1: [%.4f, %.4f]\n', principal_components(:, 1));
fprintf('PC2: [%.4f, %.4f]\n', principal_components(:, 2));
fprintf('解释方差: [%.4f, %.4f]\n', explained_variance);
fprintf('方差解释比例: [%.2f%%, %.2f%%]\n', ...
        100 * explained_variance / sum(explained_variance));

% 降维到1D
X_pca = X_centered * principal_components(:, 1);
X_reconstructed = X_pca * principal_components(:, 1)' + mean(X3);

% 绘制PCA结果
figure('Name', 'PCA分析');
subplot(1, 2, 1);
scatter(X3(:, 1), X3(:, 2), 'b.', 'MarkerSize', 8);
hold on;
% 绘制主成分方向
pc1_line = mean(X3) + 3*sqrt(explained_variance(1)) * [-principal_components(:, 1)'; principal_components(:, 1)'];
pc2_line = mean(X3) + 3*sqrt(explained_variance(2)) * [-principal_components(:, 2)'; principal_components(:, 2)'];
plot(pc1_line(:, 1), pc1_line(:, 2), 'r-', 'LineWidth', 3);
plot(pc2_line(:, 1), pc2_line(:, 2), 'g-', 'LineWidth', 2);
axis equal; grid on; title('原始数据和主成分');
legend('数据点', 'PC1', 'PC2');

subplot(1, 2, 2);
scatter(X3(:, 1), X3(:, 2), 'b.', 'MarkerSize', 8);
hold on;
scatter(X_reconstructed(:, 1), X_reconstructed(:, 2), 'r.', 'MarkerSize', 8);
axis equal; grid on; title('PCA重构');
legend('原始数据', '1D重构');

%% 实验4：图像压缩
fprintf('\n实验4：图像压缩\n');
fprintf('目标：使用SVD进行图像压缩\n\n');

% 创建简单的测试图像
[X, Y] = meshgrid(1:64, 1:64);
test_image = sin(X/10) .* cos(Y/10) + 0.5*sin(X/5) + 0.3*cos(Y/8);
test_image = (test_image - min(test_image(:))) / (max(test_image(:)) - min(test_image(:)));

fprintf('测试图像大小: %d×%d\n', size(test_image));

% SVD分解
[U_img, S_img, V_img] = svd(test_image);
s_img = diag(S_img);

% 不同压缩级别
compression_ranks = [1, 5, 10, 20, 30];
fprintf('\n图像压缩结果:\n');
fprintf('秩    PSNR(dB)    压缩比    存储节省\n');
fprintf('--    --------    ------    --------\n');

figure('Name', 'SVD图像压缩');
subplot_idx = 1;

for r = compression_ranks
    % 压缩重构
    img_compressed = U_img(:, 1:r) * S_img(1:r, 1:r) * V_img(:, 1:r)';
    
    % 计算PSNR
    mse = mean((test_image(:) - img_compressed(:)).^2);
    psnr = 10 * log10(1 / mse);
    
    % 压缩比
    original_size = numel(test_image);
    compressed_size = r * (size(test_image, 1) + size(test_image, 2) + 1);
    compression_ratio = compressed_size / original_size;
    storage_saving = (1 - compression_ratio) * 100;
    
    fprintf('%2d    %8.2f    %6.3f    %7.1f%%\n', ...
            r, psnr, compression_ratio, storage_saving);
    
    % 显示压缩结果
    if subplot_idx <= 6
        subplot(2, 3, subplot_idx);
        imagesc(img_compressed);
        colormap gray; axis equal; axis off;
        title(sprintf('秩 %d (PSNR: %.1f dB)', r, psnr));
        subplot_idx = subplot_idx + 1;
    end
end

% 显示原始图像
subplot(2, 3, 6);
imagesc(test_image);
colormap gray; axis equal; axis off;
title('原始图像');

%% 实验5：数据去噪
fprintf('\n实验5：数据去噪\n');
fprintf('目标：使用SVD进行数据去噪\n\n');

% 生成含噪声的低秩数据
true_rank = 2;
m5 = 50; n5 = 30;
U_clean = randn(m5, true_rank);
V_clean = randn(n5, true_rank);
A_clean = U_clean * V_clean';

% 添加噪声
noise_level = 0.5;
noise5 = noise_level * randn(m5, n5);
A_noisy = A_clean + noise5;

fprintf('数据矩阵: %d×%d, 真实秩: %d\n', m5, n5, true_rank);
fprintf('噪声水平: %.2f\n', noise_level);
fprintf('信噪比: %.2f dB\n', 10*log10(var(A_clean(:))/var(noise5(:))));

% SVD去噪
[U5, S5, V5] = svd(A_noisy);
s5 = diag(S5);

% 奇异值阈值选择
fprintf('前10个奇异值: [%.4f, %.4f, %.4f, %.4f, %.4f, %.4f, %.4f, %.4f, %.4f, %.4f]\n', s5(1:10));

% 硬阈值去噪
threshold = 2 * noise_level * sqrt(max(m5, n5));
fprintf('硬阈值: %.4f\n', threshold);

% 选择保留的奇异值
keep_indices = s5 > threshold;
r_denoised = sum(keep_indices);
fprintf('保留的奇异值个数: %d\n', r_denoised);

% 去噪重构
A_denoised = U5(:, keep_indices) * S5(keep_indices, keep_indices) * V5(:, keep_indices)';

% 评估去噪效果
mse_noisy = mean((A_clean(:) - A_noisy(:)).^2);
mse_denoised = mean((A_clean(:) - A_denoised(:)).^2);

fprintf('去噪效果评估:\n');
fprintf('噪声数据MSE: %.6f\n', mse_noisy);
fprintf('去噪后MSE: %.6f\n', mse_denoised);
fprintf('改善倍数: %.2f\n', mse_noisy / mse_denoised);
fprintf('PSNR改善: %.2f dB\n', 10*log10(mse_noisy / mse_denoised));

%% 实验6：推荐系统
fprintf('\n实验6：推荐系统\n');
fprintf('目标：使用SVD进行协同过滤\n\n');

% 模拟用户-物品评分矩阵
n_users = 20;
n_items = 15;
true_factors = 3;

% 生成潜在因子
user_factors = randn(n_users, true_factors);
item_factors = randn(n_items, true_factors);
ratings_complete = user_factors * item_factors' + 3;  % 基础评分3分

% 限制评分范围到1-5
ratings_complete = max(1, min(5, ratings_complete));

% 模拟缺失数据（只有30%的评分可观测）
observed_ratio = 0.3;
mask = rand(n_users, n_items) < observed_ratio;
ratings_observed = ratings_complete;
ratings_observed(~mask) = 0;  % 未观测的设为0

fprintf('用户数: %d, 物品数: %d\n', n_users, n_items);
fprintf('观测评分比例: %.1f%%\n', observed_ratio * 100);
fprintf('观测评分数: %d / %d\n', sum(mask(:)), numel(mask));

% 使用SVD进行矩阵补全（简化方法）
% 将0替换为全局平均值
global_mean = mean(ratings_observed(mask));
ratings_filled = ratings_observed;
ratings_filled(~mask) = global_mean;

fprintf('全局平均评分: %.2f\n', global_mean);

% SVD分解
[U_rec, S_rec, V_rec] = svd(ratings_filled);

% 选择合适的秩
ranks_test = [2, 3, 5, 8, 10];
fprintf('\n不同秩的预测效果:\n');
fprintf('秩    训练RMSE    测试RMSE\n');
fprintf('--    ---------    --------\n');

for r = ranks_test
    % 低秩近似
    ratings_pred = U_rec(:, 1:r) * S_rec(1:r, 1:r) * V_rec(:, 1:r)';
    ratings_pred = max(1, min(5, ratings_pred));  % 限制范围
    
    % 训练误差（观测数据）
    train_rmse = sqrt(mean((ratings_pred(mask) - ratings_complete(mask)).^2));
    
    % 测试误差（未观测数据）
    test_rmse = sqrt(mean((ratings_pred(~mask) - ratings_complete(~mask)).^2));
    
    fprintf('%2d    %9.4f    %8.4f\n', r, train_rmse, test_rmse);
end

%% 实验7：性能分析
fprintf('\n实验7：性能分析\n');
fprintf('目标：分析SVD的计算复杂度\n\n');

sizes = [50, 100, 200, 500];
fprintf('矩阵规模    SVD时间    经济型SVD    存储比较\n');
fprintf('--------    -------    ---------    --------\n');

for i = 1:length(sizes)
    m = sizes(i);
    n = min(m, 100);
    A_perf = randn(m, n);
    
    % 完整SVD
    tic;
    [U_full, S_full, V_full] = svd(A_perf);
    time_full = toc;
    storage_full = numel(U_full) + numel(S_full) + numel(V_full);
    
    % 经济型SVD
    tic;
    [U_econ, S_econ, V_econ] = svd(A_perf, 'econ');
    time_econ = toc;
    storage_econ = numel(U_econ) + numel(S_econ) + numel(V_econ);
    
    storage_ratio = storage_econ / storage_full;
    
    fprintf('%6d      %7.4f    %9.4f    %8.3f\n', ...
            m, time_full, time_econ, storage_ratio);
end

%% 实验总结
fprintf('\n=== 实验总结 ===\n');
fprintf('1. SVD提供了矩阵的最优低秩近似\n');
fprintf('2. SVD在数据压缩和去噪中表现优异\n');
fprintf('3. PCA本质上是基于SVD的降维方法\n');
fprintf('4. SVD在推荐系统中有重要应用\n');
fprintf('5. 经济型SVD在大规模问题中更实用\n');
fprintf('6. 选择合适的秩是SVD应用的关键\n\n');

end